
clear all;

%species = {'Abalone'; 'Spurpuratus'; 'Lvariegatus'; 'Lpictus';  ...
%    'Arbacia'; 'Ciona'; 'Zebrafish'; 'Bull'; ...
%    'Human-low'; 'Human-high'};
species = {'Abalone'; 'Lpictus'; 'Lvariegatus'; 'Spurpuratus';  ...
    'Arbacia'; 'Ciona'; 'Zebrafish'; 'Bull'; ...
    'Human-low'; 'Human-high'};
dates = {'200617B'; '200617B'; '200617B' ; '200617B'; ...
    '200617B'; '200617B'; '200617B' ; '200617B';...
    '200617B'; '200617B'};

% actualSpecies{1} = 'H. rufescens';
% actualSpecies{2} = 'S. purpuratus';
% actualSpecies{3} = 'L. variegatus';
% actualSpecies{4} = 'L. pictus';
% actualSpecies{5} = 'A. punctulata';
% actualSpecies{6} = 'C. intestinalis';
% actualSpecies{7} = 'D. rerio';
% actualSpecies{8} = 'B. taurus';
% actualSpecies{9} = 'H. sapiens';
% actualSpecies{10} = 'H. sapiens (viscous)';
actualSpecies{1} = 'H. rufescens';
actualSpecies{2} = 'L. pictus';
actualSpecies{3} = 'L. variegatus';
actualSpecies{4} = 'S. purpuratus';
actualSpecies{5} = 'A. punctulata';
actualSpecies{6} = 'C. intestinalis';
actualSpecies{7} = 'D. rerio';
actualSpecies{8} = 'B. taurus';
actualSpecies{9} = 'H. sapiens';
actualSpecies{10} = 'H. sapiens (viscous)';

bpWMt = [];
bpDMt = [];
grp = [];
for x = 1 : length(species)
    %newX = [7 5 4 3 2 1 6 8 9 10];
    newX = [7 3 4 5 2 1 6 8 9 10];
    
    filename = strcat(species{x},['_' dates{x} '.mat']);
    load(filename)
    
    TracktoAllSpecWMt_med_=[];
    TracktoOwnSpecWMt_med=[];
    TracktoAllSpecWMt_med=[];
    trlen=[];
    for i=1:length(procdata)
        TracktoOwnSpecWMt_med = [TracktoOwnSpecWMt_med; procdata(i).wM_track_to_median];
        TracktoAllSpecWMt_med = [TracktoAllSpecWMt_med; procdata(i).DMTracktoMedian];
        
        if ~isempty(procdata(i).wM_track_to_median)
            trlen = [trlen; size(procdata(i).A_overthresh,2)];
        end
    end
    
    %hold on; plot(trlen,selfWMt_med,'.')
    
    %need to rearrange TracktoAllSpec for some reason
    for idx = 1:10
        TracktoAllSpecWMt_med_(:,idx) = TracktoAllSpecWMt_med(:,newX(idx));
        [n(idx),h(idx)] = ranksum(TracktoAllSpecWMt_med(:,newX(idx)),TracktoAllSpecWMt_med(:,newX(x)));
    end
    
    bpWMt{x} = TracktoOwnSpecWMt_med;
    mean_bpWMt(x) = sum(TracktoOwnSpecWMt_med.*trlen)/sum(trlen);%mean(TracktoOwnSpecWMt_med);
    TEMPbpDMt{x} = TracktoAllSpecWMt_med_;
    KD_bpDMt(x,:) = sum(bsxfun(@times,TracktoAllSpecWMt_med_,trlen),1)/sum(trlen);%median(TracktoAllSpecWMt_med_,1);
    stdevKD_bpWMt(x,:) = sqrt(var(TracktoAllSpecWMt_med_,trlen/sum(trlen),1));
    grp = [grp; x*ones(length(trlen),1)];
    %ttests(x,:) = h;
    H(x,:) = h;
    P(x,:) = n;
end


%% Plot Figure S6
figure; imagesc(KD_bpDMt(1:10,1:10)); axis image; colorbar; colormap hot;
ylabel('Cell track');
xlabel('Species aggregate');
title('Weighted (by track length) Track-Species Mean Kinematic Distance')
set(gca,'xtick',[1:10],'xticklabel',actualSpecies,'ytick',[1:10],'yticklabel',actualSpecies);
xtickangle(90);

figure; imagesc(stdevKD_bpWMt(1:10,1:10)); axis image; colorbar; colormap viridis;
ylabel('Cell track');
xlabel('Species aggregate');
title('Weighted (by track length) Track-Species Mean Kinematic Distance')
set(gca,'xtick',[1:10],'xticklabel',actualSpecies,'ytick',[1:10],'yticklabel',actualSpecies);
xtickangle(90);

figure; subplot(1,2,1); imagesc(H); axis image; colorbar;
subplot(1,2,2); imagesc(P>0.05); axis image; colorbar;


%% Run donor comparison
clear all;

spec = 'Zebrafish';
date = '200617B';

filename = strcat(spec,['_' date '.mat']);
load(filename)

TracktoOwnSpecWMt_med=[];
cnt = 0;
trlen=[];
idx = [];

for i=1:length(procdata)
    cnt = cnt+1;
    TracktoOwnSpecWMt_med = [TracktoOwnSpecWMt_med; procdata(i).wM_track_to_median];
    if ~isempty(procdata(i).wM_track_to_median)
        trlen = [trlen; size(procdata(i).A_overthresh,2)];
        idx = [idx; cnt];
    end
    
    
    
end

switch spec
    case 'Spurpuratus'
        idxs{1} = 1:7; idxs{2} = 8:19; idxs{3} = 20:29;
    case 'Lvariegatus'
        idxs{1} = 1:13; idxs{2} = 14:16; idxs{3} = 17:33; idxs{4} = 34:60;
    case 'Lpictus'
        idxs{1} = 1:31; idxs{3} = 50:78; idxs{2} = 32:49; idxs{4} = 79:93;
    case 'Ciona'
        idxs{1} = 1:26; idxs{2} = 27:45;
    case 'Zebrafish'
        idxs{1} = 1:12; idxs{2} = 13:15; idxs{3} = 16;        
end

tridx = cumsum(trlen);
iii=5; jjj=5;

for n=1:length(idxs)
    tridxs{n} = tridx(idxs{n});
    trlens{n} = trlen(idxs{n});
    
    sum(TracktoOwnSpecWMt_med(idxs{n}).*trlen(idxs{n}))/sum(trlen(idxs{n}))
    sum(trlen(idxs{n}))
    
    A{n} = [];
    
    for r = 1:length(tridxs{n})
        A{n} = [A{n}, A_overthresh(:,(tridxs{n}(r)-trlens{n}(r)+1):tridxs{n}(r))];
    end
    
    
    for j =jjj:-1:1
        idx = ceil(iii*rand(size(A{n},2),1));
        
        for i = iii:-1:1
            A_ = A{n}(:,idx==i);
            
            [U_{n}, S_{n}, ~] = svd(single(A_));
            
            for k =1:2
                difft = mean(diff(U_{n}(1:3,k)));
                if difft <0
                    U_{n}(:,k) = -U_{n}(:,k);
                end
            end
        end
    end
    %sum(TracktoOwnSpecWMt_med(idxs{p}).*trlen(idxs{p}))/sum(trlen(idxs{p}))*sqrt(sum(trlen(idxs{p})))
end


sum(TracktoOwnSpecWMt_med.*trlen)/sum(trlen)

for a=1:size(idxs,2)
 for b=1:size(idxs,2)
%[a b]
KDdiffs(a,b) = abs(sum(TracktoOwnSpecWMt_med(idxs{a}).*trlen(idxs{a}))/sum(trlen(idxs{a}))-sum(TracktoOwnSpecWMt_med(idxs{b}).*trlen(idxs{b}))/sum(trlen(idxs{b})))
p(a,b) = ranksum(TracktoOwnSpecWMt_med(idxs{a}),TracktoOwnSpecWMt_med(idxs{b}))
 end
end
%%
kmat=[1 1; 1 2; 2 1; 2 2];

options=optimset('Algorithm', 'interior-point', 'TolX', 1e-8, 'TolFun', 1e-8);

for k = 1:size(kmat, 1)
    v1 = U_{1}(:, kmat(k, 1));
    v2 = U_{2}(:, kmat(k, 2));    
    c0 = [0 1 1];
    lb = [-0.1 0.5 1];
    ub = [0.1 2 1];
    options = optimoptions('particleswarm','SwarmSize',40,'UseParallel',true);
        cn = particleswarm(@(c) opt_dist(v1,-v2,c(1:2)), 2, lb(1:2), ub(1:2),options);
        cp = particleswarm(@(c) opt_dist(v1,v2,c(1:2)), 2, lb(1:2), ub(1:2),options);
    
    if opt_dist(v1, -v2, cn) < opt_dist(v1, v2, cp)
        copt{k}(1,1, :) = cn;
        nflag = -1;
        cO = cn;
    else
        copt{k}(1, 1, :) = cp;
        cO = cp;
        nflag = 1;
    end
    D(k) = opt_dist(v1, v2*nflag, cO);
    
end
        
        S(:,1) = S_{1}(1:2);
        S(:,2) = S_{2}(1:2);
        S(2,:) = S(2,:) - S(1,:);
        
        if D(1) + D(4) > D(2) + D(3)
            m1 = D(2);
            m2 = D(3);
            
            S1 = S(1,kmat(2,1)) + S(2,kmat(2,2));
            S2 = S(2,kmat(3,1)) + S(2,kmat(3,2));
            
            wM = (S1*m1 + S2*m2)/(S1+S2);
        else
            flipmat(i,j) = 1;
            m1 = D(1);
            m2 = D(4);
            
            S1 = S(1,kmat(1,1)) + S(2,kmat(1,2));
            S2 = S(2,kmat(4,1)) + S(2,kmat(4,2));
            
            wM = (S1*m1 + S2*m2)/(S1+S2);
            
        end
        